#R version 4.2.1 (2022-06-23)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19044)

library(dplyr) #version 1.0.9
library(readr) #version 2.1.2
library(stringr) #version 1.4.0
library(zoo) #version 1.8-10
library(sandwich) #version 3.0-2
library(lmtest) #version 0.9-40
library(prediction) #version 0.3.14
library(ggpubr) #version 0.4.0
library(estimatr) #version 1.0.0
library(margins) #version 0.3.26
library(broom) #version 1.0.0
library(scales) #version 1.2.0
library(stargazer) #version 5.2.3

`%notin%` <- function(x,y) !(x %in% y) 

a1 <- data.frame(xpos = c(-Inf,-Inf,Inf,Inf), ypos =  c(-Inf, Inf,-Inf,Inf), annotateText = c("GWF"), hjustvar = -0.25, vjustvar = -0.5)
a2 <- data.frame(xpos = c(-Inf,-Inf,Inf,Inf), ypos =  c(-Inf, Inf,-Inf,Inf), annotateText = c("DD"), hjustvar = -0.25, vjustvar = -0.5)
a3 <- data.frame(xpos = c(-Inf,-Inf,Inf,Inf), ypos =  c(-Inf, Inf,-Inf,Inf), annotateText = c("PAR"), hjustvar = -0.25, vjustvar = -0.5)
a4 <- data.frame(xpos = c(-Inf,-Inf,Inf,Inf), ypos =  c(-Inf, Inf,-Inf,Inf), annotateText = c("WTH"), hjustvar = -0.25, vjustvar = -0.5)

df <- read_csv("data/gwf2010.csv") %>% 
  select(gwf_regimetype, year) %>% 
  filter(!is.na(gwf_regimetype), year >=1942) %>% 
  mutate(gwf_regimetype = ifelse(gwf_regimetype %in% c("indirect military", "military", "military-personal", "party-military", "party-military-personal"), "Military", "Civilian")) %>%
  filter(gwf_regimetype %in% c("Civilian", "Military")) %>% 
  group_by(year) %>% 
  count(gwf_regimetype) %>% 
  mutate(percent = n/sum(n)) 

gwf <- ggplot(df) +
  geom_line(aes(x = year, y = percent, color = gwf_regimetype)) +
  geom_point(aes(x = year, y = percent, shape = gwf_regimetype)) +
  geom_text(data=a1,aes(x=xpos,y=ypos,hjust=hjustvar,vjust=vjustvar,label=annotateText)) +
  theme_bw() +
  xlab("") +
  ylab("Regime Frequency") +
  scale_color_grey() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits=c(0,1)) +
  guides(color = guide_legend(title = "Regime"), shape = guide_legend(title = "Regime")) +
  theme(legend.position="none") 

#dd case
df <- read_csv("data/DD.csv") %>% 
  filter(democracy == 0) %>% 
  dplyr::select(countryname, year, emil) %>% 
  filter(!is.na(emil)) %>% 
  group_by(year) %>% 
  count(emil) %>% 
  mutate(p = n/sum(n)) %>% 
  na.omit() %>% 
  mutate(emil = recode(emil, `1` = "military", `0` = "civilian")) 

dd <- ggplot(df) +
  geom_line(aes(x = year, y = p, color = emil)) +
  geom_point(aes(x = year, y = p, shape = emil)) +
  geom_text(data=a2,aes(x=xpos,y=ypos,hjust=hjustvar,vjust=vjustvar,label=annotateText)) +
  theme_bw() +
  xlab("") +
  ylab("") +
  scale_color_grey() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits=c(0,1)) +
  guides(color = guide_legend(title = "Regime"), shape = guide_legend(title = "Regime")) +
  theme(legend.position="none") 

#svolik
df <- read_csv("data/Svolik clean.csv") %>% 
  mutate(mil = recode(par_military, "corporate" = "Military", "indirect" = "Military", "personal" = "Military", "civilian" = "Civilian")) %>% 
  dplyr::select(countryname, year, mil) %>% 
  filter(mil %in% c("Civilian", "Military"), year %in% 1946:2008) %>% 
  group_by(year) %>% 
  count(mil) %>% 
  mutate(p = n/sum(n)) 

par <- ggplot(df) +
  geom_line(aes(x = year, y = p, color = mil)) +
  geom_point(aes(x = year, y = p, shape = mil)) +
  geom_text(data=a3,aes(x=xpos,y=ypos,hjust=hjustvar,vjust=vjustvar,label=annotateText)) +
  theme_bw() +
  xlab("Year") +
  ylab("Regime Frequency") +
  scale_color_grey() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits=c(0,1)) +
  guides(color = guide_legend(title = "Regime"), shape = guide_legend(title = "Regime")) +
  theme(legend.position = "bottom") 

#wth
df <- read_csv("data/ARD.csv") %>%
  filter(regimeny %notin% c("democracy", "civil war", "occupation", "other", "rebel regime", "transitional"), !is.na(regimeny)) %>% 
  mutate(mil = ifelse(str_detect(regimeny, "military"), "Military", "Civilian")) %>% 
  dplyr::select(year, mil) %>% 
  filter(mil %in% c("Civilian", "Military")) %>% 
  group_by(year) %>% 
  count(mil) %>% 
  mutate(p = n/sum(n)) 

wth <- ggplot(df) +
  geom_line(aes(x = year, y = p, color = mil)) +
  geom_point(aes(x = year, y = p, shape = mil)) + 
  geom_text(data=a4,aes(x=xpos,y=ypos,hjust=hjustvar,vjust=vjustvar,label=annotateText)) +
  theme_bw() +
  xlab("Year") +
  ylab("") +
  scale_color_grey() +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits=c(0,1)) +
  guides(color = guide_legend(title = "Regime"), shape = guide_legend(title = "Regime")) +
  theme(legend.position = "bottom") 

ggarrange(gwf, dd, par, wth, ncol =2 , nrow = 2) 

#save to file
ggsave(filename = "figures/fig1.jpg", dpi = 500, width = 5, height = 5)
